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ABSTRACT 

A numerical study of an algorithm proposed by Gusein Guseinov, which determines approx- 
imations to the optimal solution of problems of calculus of variations using two discretiza- 
tions and correspondent Euler- Lagrange equations, is investigated. The results we obtain 
to discretizations of the brachistochrone problem and Mania example with Lavrentiev's phe- 
nomenon are compared with the solutions found by other methods and solvers. We conclude 
that Guseinov's method presents better solutions in most of the cases studied. 

Keywords: calculus of variations, Euler-Lagrange equations, discretization, solvers, brachis- 
tochrone problem, Lavrentiev phenomenon. Mania example. 
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1 Introduction 

In a problem of the calculus of variations, functions that extremize a given functional are 
sought. Several methods to determine approximated solutions for such problems are known 
in the literature. Here we do a comparison study of different approaches. More precisely, 
we investigate a method proposed by Guseinov in [7\ that determines approximated solutions 
to problems of the calculus of variations by discretization. The results we obtain applying 
Guseinov's method to well known problems are then compared to the results found by other 
methods. 



2 Calculus of variations 

Given real numbers Xa, Xb, a, (3, Xa < Xb, and a function / : [xa,^;,] x M x M — > R called the 
Lagrangian, f £ C^, the fundamental problem of the calculus of variations aims to determine 

*This is a preprint of a paper whose final and definite form appeared in Int. J. Math. Stat. 9 (2011), no. Sll, 
26-41. 
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{x,y{x),y'{x)) = — ( — {x,y{x),y'{x))] . (2.3) 



y E C^{[xa, Xb] ) such that function y minimizes the functional 

J[y{-)] = / f{x,yix),y'{x))dx (2.1) 
under the given boundary conditions 

y{xa) = a, y{xb)=P. (2.2) 

A necessary condition for function y to solve problem (|2.ip - (|2.2|) is given by the Euler-Lagrange 
equation 

d fdf 

dy 

A solution y of ()2.3p is said to be an extremal, meaning that y either minimizes, maximizes, or 
acts like a "saddle function" of 

Example 2.1 (the brachistochrone problem). Given two points A = (xa,ya) and B = {xf),yb) 
in the same vertical plane, we want to determine the curve described by a particle that, with 
zero initial speed and under action of gravity, connects the points in a minimum time, ignoring 
friction: 

. 1 r" /i+y'(x)2 , ^ 

min^ / J \ \ dx, (2.4 

vi-) A. ^ Va- y{x) 

where g is the gravitational constant. 

The solution to problem (|2.4|) is an arc of cycloid: the curve described by the revolution of a 
point belonging to a circle of radius r that rolls without sliding on the straight line y = ya- 
Such curve is described by the parametrization 



X = Xa + r {9 — sin(6')) 
y = Va - r {1 - cos{9)) 



%<e< 01. 



Since we know {xa,ya) and {xb,yb), it is possible to determine r, Oq and 9i. Using the point 
{xa-iVa) and considering r / 0, = is obtained. The time it takes the particle to travel the 
curve is given by 



1 , , ILW' 

^ ^ \ dx I 



1 I" 1 1 r\ _\vii)j_dx 



For instance, the brachistochrone problem ()2.4p with boundary conditions 

Xa = 0, ya = W, Xb = W, yb = (2.6) 
has the following function as optimal solution: 

f j; = 5.72917061 - 5.729170 sin(^) 

< ^ , < 6l< 2.412011. 

[ y = 4.270830 + 5.729170 cos(6') 

Considering (j2.5p . the optimal value, disregarding the constant is 



r-2.412011 



^ ^ I 32.823393 sin(6>) ^ 



, 5.729170-5.729170 cos(6») , , ,^ 

\ ^ (5.729170 - 5.729170 cos(0)) dO = 8.164699. 

\ 5.729170 - 5.729170 cos(6') ^ ^ " 
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Figure 1: Solution to the brachistochrone problem with A = (0, 10) and B = (10, 0). 
Example 2.2 (Mania's example). The problem consists of finding the function y that minimizes 

I[yi-)]= (\y{xf-xf{y'{x)fdx (2.7) 

JO 

under the boundary conditions y(0) = and y(l) = 1. Mania example is a problem of the cal- 
culus of variations that exhibits the so called Lavrentiev phenomenon [21 0] . Problems with the 
Lavrentiev phenomenon have different solutions in the space of absolutely continuous functions 
and in or the space of Lipschitzian functions: the infimum of the functional in the set of 
absolutely continuous functions is strictly less than the infimum of the same functional in the ad- 
missible set C^. This gap makes it harder (or even impossible) to determine the minimum of the 
variational functional computationally. The optimal solution for the Mania problem is the func- 
tion y{x) = xA. Indeed, for all x G [0, 1] one has L {x,y{x),y' (x)) = [y{x)^ — x)^ (y'(x))^ > 0, 
and thus I[y] > 0. Since L {x,y{x),y' (x)) = (^y{x)^ — x)'^ (y' (x))^ = for all x G [0,1], then 
= is the optimum value for (|2.7|) . Although the optimal solution for this problem is 
simple, it is an open question how to find the optimal value zero for I by means of numerical 
methods. 

3 Discrete-time calculus of variations 

A way to find approximate solutions to problems (|2.ip - ()2.2p consists in subdividing the problem 
into "smaller" problems easier to solve. With this in mind, the interval is considered as 

a union of n intervals: 

[Xa,Xb] = [Xa,Xi] U [xi,X2] U • • • U [x„_i,Xfc] . 

Considering all the subintervals with the same amplitude, 

n 

[Xa,Xb\ = [J[xi-i,Xi\ , (3.1) 
1=1 
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xq = Xa, Xn = Xb, Xi = xq + zAx, z = 1, 2, . . . , n — 1, and Ax = ^>> ^ y/ approximated by 
the slope of the hne defined by the two extreme values of each interval: 



y'{xi) 



y{xi) - y{xi-i) _yi- yi-i _ ^yi-i 



Ax Ax Ax 

3.1 The standard discretization 

There are several ways to formulate the problem using the intervals (|3.ip . The most common 
approach was presented by Euler himself, and we call it the "standard" discretization, also 
known as Euler's method of finite differences. In this standard approach one searches y^, i = 
l,2,...,n — 1, that minimize (or maximize) the finite sum 



i=l ^ ' 



The Euler-Lagrange equations for this formulation are 



^ rtl' {Xi,yi, 



df f Ayi^i\ _ ^dy' \xi,yi, 
t;~ 1 Xi,yi, 



dy \ Ax ) Ax 

i = 1, 2, . . . , n — 1. 

Example 3.1 (brachistochrone problem). Using the standard discretization, the brachistochrone 
problem is approximated as follows: 



mm 

1=1 



2 

1 + > A.,. 



Va - Vi 

subject to y{xa) = ya , y{xb) = yb- 



Example 3.2 (Mania's example). The standard discretization of the Mania example is given by 



1=1 

subject to y(0) = , y(l) = 1. 



3.2 Guseinov's discretization 

Guseinov proposes in [2| a different discrete approach to find solutions to problems of the calculus 
of variations. Considering a finite set of n + 1 real numbers from the integration interval, 

X = {xo,xi, ...,Xn}, 

Xa = Xq < xi < • • • < Xn = Xb, another set of n + 1 real numbers is found: 

{yo = y{xo) = a,yi = y(xi), ...,yn = y{xn) = P}- 

Then, it is possible to define a polygonal line that connects all the points 

(xo,a), (xj,yi), (xrt,/3) , i = l,...,n - 1. 
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Guseinov's discretization is different in the way the functional is determined. The discrete 
functional used by Guseinov [7] is obtained considering the meaning of the problem and its 
implications when we discretize it. For instance, if we search a curve that extremizes a func- 
tional, such line is approximated using several straight line segments. If an area is searched, 
then we use trapezoids to define the discrete functional; when we search a revolution solid we 
approximate it with frustums of cones. The discrete Lagrangian of Guseinov is defined by 

L: XxXxMxMxM — > M 

{s,t,u,v,w) I—)- L{s,t,u,v,w), 

where ^ , ^ , ^ G C'' [M] and the variables s and t concern the domain of the functional (|2.ip , 
u = y{s), V = y{t), and w = jE^- The finite sum to minimize (or maximize) is given by 

J[y] ^Yl^ (^Xi-i,Xi, yi-i,yi, Axi_i, (3.2) 

where Ayj_i = yi — yi-i and Axj-i = Xi — Xi-i. The function y to be found verify the boundary 
conditions and should minimize (or maximize) J^. 

Example 3.3 (brachistochrone problem). To discretize the brachistochrone problem a la Gu- 
seinov [7J , we consider that the solution curve will be approximated by a sequence of n straight 
line segments that connect the points (j;j_i, yj_i) and {xi,yi), i = 1, . . . ,n, Xi > > yi. 

The amount of time spent by a particle to travel each of these line segments (neglecting friction) 
should be calculated. Regarding a, as the acceleration of the particle and 6i as the angle formed 
with the x's axis of the segment i, a, = gsinOi. By simple trigonometry, 

yi~i - yi ^y-i-i 



sml 



{xi - Xi^if + (yi_i - yif \J (Axj_i)^ + (-Ayi_i)^ 



Thus, ai = , ^^f-^ Because a^ = ^,^ = , ^^f then 

g^yi-i , , 
Vi = — =t + c. 

(Axi„i)2 + {/^yi-if 

Besides, if t = 0, then c = vq. In each line segment the counting of time is resumed (we count 
the time needed to run each line segment separately, adding all the times at the end) . This way, 
in each line segment, ^initial ~ ^' Thus, c = Vj-i, = 0, and 

g^yi-i 



{Axi^if + {/\yi-if 



--t + u 



i-i- 



Moreover, v = and therefore 



vdt = I ds, 
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where ti is the time the particle takes to run the hne segment i: 



ti 



t + Vi-i at 



V V(Ax,_i)^ + (Ay,_i)^ 

/ tdt + Vi-i 
Jo Jo 



^(Ax,_i)2 + (Ay,_if 



dt 



g^Vi-i ,2 , . 



2^(Axi_i)" + (Ay,_i)^ 



To determine the hne integral ds, a parametrization of the line segment which connects the 
points {xi-i,yi-i) and {xi,yi) is used: 



t G Xj], and 



x{t) = t, 

Jxi-i 



2 

dt 



'i+i^ydt 



Then, 



Axi_i 



Using the quadratic formula, 



VAxj_i 



2 



^ gAj/i^l 2 . , A /l _L Ayi_i 

^ , =ti - Vi-iti + Axi-i\ 1 + [ =0. 



= — ^A^~]^ — ( ""^-^ ^ ^/''^-i " '^^^y^-' I ■ ^^-^^ 



Since Ayj_i < and A.Tj_i > 0, then gAy i ' ^ Moreover, since the amount of 

time must be non-negative, the condition 

Vi-i ± ^ vf_i - 2gAyi^i < 
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must be true and 



-'i-i 



=^ Vi-i > yvf^i - '^g^Vi-i (since Vi-i > 0) 
Vi-i - - 2gAyi^i > . 

Since the expression should depend only on the points used to discretize the curve, the initial 
velocity of the particle in each line segment i (fi-i) shall be determined using these values. 
Applying the energy conservation law (for the importance of conservation laws in the calculus 
of variations we refer the reader to O [6] ) , the amount of mechanical energy in the beginning 
of the first straight line segment (Em^) is the same as the amount of energy at the beginning 
of each of the other line segments. Considering the line segment i, = Efm- At any given 
moment, the mechanical energy of the particle is the sum of its potential and kinetic energies. 
Considering g as the gravity acceleration, m the mass of the particle, h the height at which the 
particle is, and v the velocity of the particle, then 

Ep = gmh and E^ = ^ . 
Therefore, in the line segment i (with extremes at (2;j_i,yj_i) and {xi,yi)), 

Ema = Errn <^Ep^ + Ec^ = Ep^ + Em^ 

■^gmya + = gmyi^i -\ ^ — 

■^gva + y = avi-i + 

Since the particle starts still, i.e., vq = 0, 



gva = gvi-i + -hr ^ "^-i = ^V'^g {va-Vi-i 



Besides, Vi-i > for any i = 1,2, ...,n, so = ^J2g {ya — yi~i)- Replacing this last 

expression for Vi-i into (|3.3p . the amount of time that a particle takes to run the line segment 
i is determined: 



c/Ay~^ \ ^'^^ ~ ^'^^^ ~ iVa - Vi-i) - 2g {yt - yi-i)j 

{\/2g [ya-yi-i) - \/2g [ya-yi]^ 



Axi_n/1 + 



g^Vi-i 

- (\/ya - Vi-i - y/Va- Vi) Axi_i 



g^y^ 



V2^l+ (s^Ty) iVVa - Vi-l - VVa - Vi) ■ WVa " Vi-l + VVa " Vi) AXj-i 



1 + 



g ^Vi-l WVa - Vi-l + y/Va - Vi 



^/g^Vi-l iVVa - Vi-l + VVa - Vi) 

{ya - yi-1 -ya+ Vi) Axi_i 



A \ 2 
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1 _L ( ^y^-^ 



2 

-Vi-l + Vi) Axi_i 



1 + I ^ 



Adding the times of the n hne segments we obtain 



To find an approximated solution for the brachistochrone problem, the y^'s are searched so that 



2 

1+ 



Ay»-i 

EV yAXi_ly 
— Axj-i. 



The problem may, then, be formulated in the discrete (Guseinov) form as 



1 + 



mm > — A.T; 

subject to y{xa) = Va , y{xb) = Vb- 

Remark 3.1. Since the Mania example is a theoretical problem, without no a priori physical 
meaning, it is not clear how to discretize it using Guseinov's approach. 

3.3 Euler-Lagrange equation 

The next result presents the Euler-Lagrange equation for the discrete Guseinov formulation of 
the problem of the calculus of variations. 

Theorem 3.1 (Guseinov's discrete Euler-Lagrange equation [7J). // the variational functional 
of the calculus of variations (|3.2|) has a local extreme in y, then y satisfies the Euler-Lagrange 
equation 

dL f Ay A Axi , dL f Ay,_i 
Xi,Xi+i,yi,yi+i, - — h Xi-i,Xi,yi_i,yi 



(3.4) 



ou \ Axi J Axi^i ov \ Axi-i 



Axi^l 

forie {1,2,... ,n- 1}. 

Proof See [7]. □ 

Remark 3.2. The Euler-Lagrange equation (|3.4p is to be complemented with the given boundary 
conditions y{xQ) = a and y{xn) = /?• 
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Using Theorem l3.1l fand the given boundary conditions), Guseinov [7] developed a (new) method 
that, discretizing the interval, determines approximated solutions for problems of the calculus 
of variations. We note that Theorem 13.11 does not imply the use of a "Guseinov discretization". 
There is no reference to the kind of discrete functional L. Functional L must, only, verify the 
conditions ^ £ C°[M]. In other words, the method may be applied using the "stan- 

dard discretization" as long as the standard discrete functional verifies such conditions. The 
approximations to solutions found by application of Theorem 13 . 1 1 are candidates to extremize the 
functional (they are critical functions), which means that the found functions may maximize, 
minimize or neither maximize nor minimize the functional {saddle function). 

3.4 Guseinov's algorithm 

Based on Theorem 13.11 the next algorithm determines approximated solutions for problems of 
the calculus of variations. 

Input: a continuous Lagrangian /, a discrete Lagrangian L, boundary conditions y{xa) = a 
and y{xh) = /3, and the number n of intervals that divide the integration domain. 

Output: y{xi), i = 0, 1, . . . ,n, J[y] given by (j3.2p . and / {x,y{x),y' (x)) dx, where y is the 
piecewise linear function defined by the points y{xi) found by the algorithm. 

Algorithm 

1. Read input data: L{s,t,u,v,w), to ^ Xa, tn ^ Xf,, a ya, (3 ^ yt, and n. 

2. Determine the partial derivatives of L: Lu ^ ^ {s,t,u,v,w), Lv ^ ^ {s,t,u,v,w), 

^ ^ {s,t,u,v,w). 

3. Determine the set of n + 1 equidistant x: 

• Ax = SillSs^; 

• X = {xq = Xa, xi = xq + Ax, . . . ,Xi = xo + iAx, . . . , x„ = Xb}- 

4. Solve the discrete (nonlinear) Euler-Lagrange equations system (|3.4p : 

• determine the system of equations F{Y) = such that 

— for i = 1, . . . , n — 1, 

Fi = Lu (^Xi,Xi+i,yi,yi+i, ^'+^~^' ^ + Lv (^Xi.i,Xi,yi-i,yi, 

Ax 

— the boundary conditions Fq = yo — a, Fn = y-n — (3 hold; 

• solve the system. 



5. Compute J[y{-)] = '^L fxi^i,Xi,y{xi_i),y{xi), M Ax. 

1=1 ^ ^ ^ 
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6. Determine the piecewise linear function y. 



7. Compute the value of J^^ f{x,y{x),y'{x))dx. 

8. Draw the graphic of function y. 

3.5 Implementation 

The algorithm of ^3.41 was implemented using Mathematica 6. Considering the problem and 
the possible amount of points chosen from the domain of the integral, the nonlinear system of 
equations is not, in general, easily solved. So, in general, the Mathematica' s functions used to 
solve the system of equations may not find solutions or may find complex numbers. Neither of 
these results is acceptable. In our implementation, a function that numerically solves nonlinear 
programming problems is used. In this way an approximation to the solution of the system of 
equations is found. Given the nonlinear equations system 

F{i) = 0, i = 1,... ,n- 1 

< yo - a = 
. yn-/3 = 

we consider the following nonlinear programming problem: 

n-l 

min ^(F(0)2 + (yo-a)2 + (yn-/3)' 

^=1 (3.5) 
subject to yo = a, yn = /3 , 

To find a numerical approximation to the solution of the nonlinear programming problem (|3.5|) , 
the function NMinimize of Mathematica is used. It is possible to choose the method used to 
find solutions: Random Search, Nelder Mead, Differential Evolution, or Simulated Annealing. 
These methods present solutions that may or may not be the same. Besides, one of the solutions 
may be the best considering the nonlinear problem (solving the nonlinear equations system) but 
not the best regarding the initial problem. So, from the presented solutions by these methods, 
the one that minimizes the discrete problem of the calculus of variations is the one chosen. 
Although some methods end with success, they may present solutions that make no sense (such 
as yi close to plus or minus infinity). Thus, since the nonlinear programming problem may have 
several restrictions, a new set of restrictions was added, trying to improve the solutions. This 
way, two set of restrictions were tested: 

• "restrictions 1" : boundary conditions from the problem of the calculus of variations; 

• "restrictions 2" : boundary conditions from the problem of the calculus of variations and 

< max{|a|, |/3|}, i = l,...,n-l. (3.6) 
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3.6 An optimal control solver 

OC - Optimal Control solver - is a solver supplied in [8]. This solver determines approximated 
solutions for optimal control problems in Mayer formulation. The solutions are found by dis- 
cretizing the interval [x^jXt] in n subintervals and using Mathematical Programming methods. 
It is possible to choose the method used to optimize {Conjugated Gradients - CG, Newton 
Method - NM, Univariate Search - US, Direct Search ~ DS 1 and DS 2, and Random Search 
- RS), the method used to solve the differential equations {Euler or Runge-Kutta), the initial 
solution, maximum number of iterations, and the number of intervals of the discretization. To 
use this solver, the problem of the calculus of variations must be formulated as an optimal 
control problem: 

min / L {x,y{x), z{x)) dx 

subject to y'{x) = z{x) 

y{xa) = a, y{xb) = 13. 

This optimal control problem in Mayer form is: 

min u{xb) 
subject to y'{x) 
u'{x) 
u{xa) 

y{xa) 

3.7 An evolutionary algorithm 

A simplification of the (^, A)— ES Algorithm (an evolutionary algorithm) presented in ^j, uses 
evolutionary strategies combined with optimal control to find approximated solutions to prob- 
lems of the calculus of variations. The algorithm keeps seeking solutions to the problem, eval- 
uating each of them. The solutions closer to the target set are used to find new solutions 
(supposedly better). This process ends after a certain number of iterations (see [3] for details). 

4 Results and comparisons 

The result to be compared is the value of the functional integral (|2.ip along the approximated 
solutions. 

Example 4.1 (brachistochrone problem). Table [T] presents the value of the integral (|2.5p using as 
integrand the piecewise linear functions defined with the approximated solutions of (|2.4p found 
by the Guseinov algorithm of ^3.41 and the piecewise linear function defined with the optimal 
solution (PLFOpt), whose value is used as reference. 

The best result found by our implementation of the algorithm used Guseinov's discretization, 
20 intervals for discretization, and the Random Search method (Figure [2]). The integral value 
is 8.189. Using the same options except the type of discretization - standard discretization - 
the value of the integral is 8.751, which shows that the approximation is clearly worse. 



= z[x) 

= L{x,y{x),z{x)) 
= 

= a, y(xb) = f3. 
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Table 1: Values of J^'' f{x,y{x),y'{x))dx for problem ()2.4p subject to Xa = 0, x^ = 10, Ua = 10, 
and Uh = 0. 



No. 




Guseinov discretization 


Standard discretization 


Integral 


of 


Method 


Restrictions 1 


Restrictions 2 


Restrictions 1 


Restrictions 2 


PLFOpt 


Ints 






(fTHD and (13:6]) 








3 


RS 


8.353 


8.353 


8.418 


8.418 


8.369 




NM 


8.353 


8.353 


8.418 


8.418 






DE 


* 


* 


* 


8.418 






SA 


8.353 


8.353 


8.418 


8.418 




5 


RS 


8.271 


8.271 


8.464 


8.464 


8.281 




NM 


8.271 


8.271 


8.464 


8.464 






DE 


* 


* 


* 


* 






SA 


8.271 


** 


8.464 


8.464 




8 


RS 


8.229 


8.229 


8.566 


8.566 


8.235 




NM 


8.229 


9.037 


8.566 


9.906 






DE 


* 


* 


* 


* 






SA 


*** 


8.229 


8.566 


12.288 




10 


RS 


*** 


12.716 


8.617 


12.682 


8.220 




NM 


8.216 


12.716 


8.617 


8.617 






DE 


* 


* 


* 


* 






SA 


*** 


12.716 


8.617 


15.370 




15 


RS 


*** 


** 


8.702 


12.749 


8.201 




NM 


*** 


12.764 


8.702 


8.702 






DE 


* 


* 


* 


* 






SA 


* 


* 


8.702 


* 




20 


RS 


8.189 


** 


8.751 


12.782 


8.191 




NM 


*** 


15.972 


8.751 


21.053 






DE 


* 


25.115 


* 


* 






SA 


*** 


* 


* 


* 





'*': method used to solve the nonlinear programming problem doesn't end successfully (e.g., 
the method exceeded the number of iterations); 
'**': the value found is complex; 

'***': the value found is too high and the solution makes no sense; 

Methods: RS — Random Search, NM — N elder Mead, DE — Differential Evolution, 

SA — Simulated Annealing. 

Considering the mentioned algorithms, solver, and the solutions determined by them, it is 
possible to verify which one presents the best solution regarding the value of the integral - see 
Table E 
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2 4 6 S 10 

Figure 2: Approximation found applying Guseinov's discretization with n = 20 intervals and 
the Random Search method; and the exact solution. 



Table 2: Values of the integral along approximated solutions obtained by different methods. 



Solution 


Value of Integral 


Notes 


Optimal 


8.16470 




Points from the optimal line 


8.19139 


21 points with 
(successive) equidistant abscissae 


Guscinov algorithm 


8.189344 


20 intervals 
Guseinov discretization 
Method: Random Search 


Evolutionary algorithm 


8.19365 


20 intervals 


OC 


8.336 


20 intervals 
Method: Conjugated Gradients 
Method: Runge-Kutta 2 
Piecewise Linear 



Example 4.2 (Mania's example). The results obtained with our implementation of Guseinov's 
algorithm ( ^3.4p using the standard discretization with two different sets of restrictions in the 
nonlinear programming problem and different methods are presented on Table El 

The solution of the best "numerical" result obtained using our implementation in Mathematica 
applied 8 discretization intervals and the method of Simulated Annealing (Figure [3]). The result 
is the value 0.328 for the integral ()2.7p . 

Approximations for the solution found with n = 10, n = 15 or n = 20 discretization intervals 
seem to be closer to the optimal curve than de previous one (Figured]). However, the results 
with n = 10, n = 15 or n = 20 are worse than with n = 8 (for instance, the value of the 
functional obtained with 15 intervals is 1.156). 

Moreover, taking a close look at the last column of the Table [3l something is apparently wrong 
because the value of the integral is increasing as the number of discretization intervals increases, 
instead of becoming closer to zero (the optimal value). Although this looks contradictory, it 
is not. In fact this is a consequence of the Lavrentiev phenomenon exhibited by this problem. 
The next result proposed and proved in [1] explains the fact. 
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Table 3: Values of J^*" L{x,y{x)^y' {x))dx for the functional (|2.7p {xa = 0, x;, = 1, 2/a = and 

yb = i)- " 



Number of 


Method 


Standard discretization 


Integral 


Intervals 




Restrictions 1 


Restrictions 2 


PLFOpt 


3 


RS 


92.314 


92.314 


0.229 




NM 


92.314 


92.314 






DE 


92.314 


92.314 






SA 


92.314 


92.314 




5 


RS 


1488.100 


1488.100 


0.381 




NM 


1488.090 


0.994 






DE 


1488.100 


1488.100 






SA 


1488.100 


1488.100 




8 


RS 


17549.400 


14.768 


0.610 




NM 


539.515 


537.198 






DE 


64.207 


17549.400 






SA 


17549.400 


0.328 




10 


RS 


55619.100 


0.490 


0.762 




NM 


55619.100 


20.364 






DE 


2390.020 


55619.000 






SA 


55619.100 


3.750 




15 


RS 


443732.000 


1.156 


1.143 




NM 


443732.000 


12416.200 






DE 


0.474 


443732.000 






SA 


443732.000 


66.265 




20 


RS 


1915810.000 


127.346 


1.524 




NM 


1915810.000 


10.174 






DE 


2235.890 


32784.700 






SA 


0.956 


7.473 






Figure 3: Approximation obtained by applying the standard discretization with "Restrictions 2", n = 8, and the 
Simulated Annealing method; together with the optimal absolute continuous curve. 
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0.2 OA 0.6 O.S 1.0 

Figure 4l Approximation obtained by applying the standard discretization with "Restrictions 2", n = 15, and the 
Random Search method. 

Theorem 4.1. For any sequence of Lipschitz trajectories {yn}„ such that y„ tends to y{x) = xa 
as n tends to oo, for almost all x € [0, 1], then T[yn] defined in (|2.7p tends to oo. 

Considering the peculiarity of the Mania example, it is very difficult to compare objectively the 
results determined by other methods and algorithms. Unlike the brachistochrone problem, there 
is apparently no relation between the solution that is graphically better and the solution that 
is numerically better. However, regarding the graphics, the solutions found by the Guseinov 
algorithm are good. Increasing the number of intervals used in the discretization, the approxi- 
mations improve considering the graphical representation but the integral value becomes worse. 
This fact is also seen in the graphics of the piecewise linear functions defined using points of the 
optimal solution and in the approximations found by the OC solver. The best approximation 
found by the OC solver, considering the value of the integral, used the methods of Conjugated 
Cradients, Runge-Kutta, and Piecewise Constant, with 0.0326998 as the value of the integral. 
However the graphic of the approximation found using the methods Univariate Search, Runge- 
Kutta, and Piecewise Linear, shows that this approximation is closer to the optimal solution, 
although its integral value is 1.45193. 

In [1] the Truncation Method is proposed. This method determines an upper limit for the 
integral minimum using an auxiliary functional whose integration domain is in the original 
integration domain and does not include the points that "create" the Lavrentiev phenomenon. 
This new integration domain is also divided in n intervals. The numerical results presented in 
the paper [1] don't include the integral value used in our work as the key comparison element. 

5 Conclusion 

The algorithm proposed by Guseinov in [7] presents very good solutions to "regular" problems 
of the calculus of variations. Moreover, the usage of Guseinov discretization improves the 
results. The method works worse when applied to problems with the Lavrentiev gap, like the 
Mania example, or to problems of optimal control. Since the method involves the resolution of 
nonlinear systems of equations, if the numerical methods used in the solvers are not working 
well the algorithm does not find solutions as good as it could. The solutions obtained for the 
Mania example are not worse than the solutions found by other methods and solvers. Our 
implementation of the algorithm in Mathematica (version 6) is very easy to use. 
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